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ABSTRACT 

A  new  adaptive  finite-difference  scheme  for  scalar  hyperbolic  conserva¬ 
tion  laws  is  introduced.  A  key  aspect  of  the  method  is  a  new  automatic  mesh 
selection  algorithm  for  problems  with  shocks.  We  show  that  the  scheme  is 
L^-stable  in  the  sense  of  Kuznetsov,  and  that  it  generates  convergent 
approximations  for  linear  problems.  Numerical  evidence  is  presented  that 

indicates  that  if  an  error  of  size  c  is  required,  our  scheme  takes  at  most 
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0(e  )  operations.  Standard  monotone  difference  schemes  can  take  up  to 

-4 

0(e  )  calculations  for  the  same  problems. 
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SIGNIFICANCE  AND  EXPLANATION 


Certain  problems  in  gas  dynamics,  oil  reservoir  simulation  and  other 
fields  can  be  modeled  by  hyperbolic  conservation  laws,  a  class  of  partial 
differential  equations.  The  solutions  of  such  problems  are  typically  made 
of  smooth  surfaces  separated  by  discontinuities,  or  shocks.  Usually, 
less  information  is  needed  to  specify  the  solution  in  the  smooth  regions 
than  in  the  shock  regions. 

1  c  .  r 

In  this  paper  we  introducers  stable  finite-difference  scheme  for  con¬ 
servation  laws  that  incorporates  a  time-varying,  nonuniform  computational 
mesh.  At  any  given  time,  ear  mesh  selection  algorithm  chooses  a  mesh  based 
on  the  approximation  calculated  up  to  the  time.  The  algorithm  uses  know¬ 
ledge  of  a  solution's  structure  to  reduce  the  number  of  meshpoints  in  the 
regions  where  the  solution  is  smooth.  This  reduces  the  method's  computa- 

■'j  _  /•.  ■  .■ 

tional  complexity  while  maintaining  full  accuracy.  Me  prove^that  mnr 
method  is  stable  for  the  complete  nonlinear  problem,  and  that  it  converges 
for  linear  problems.  -We  qiwe  examples  where  ear  method  is  asymptotically 
faster  than  previous  ones. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
stawsary  lies  with  MRC,  and  not  with  the  author  of  this  report. 


A  Stable  Adaptive  Numerical  Scheme  for  Hyperbolic 
Conservation  Laws 
* 

Bradley  J.  Lucier 

1.  Introduction. 

Our  focus  in  this  paper  is  the  efficient  solution  of  the  hyperbolic  conserva¬ 
tion  law 

xeR  *>0,  (C) 

tt(x.O)  =  u0(x). 

We  use  an  adaptive  finite-difference  scheme  that  takes  advantage  of  the  struc¬ 
ture  of  the  solution  of  (C)  to  reduce  its  computational  complexity.  We  consider 
the  scalar  case  u(x,t)cR  For  the  general  nonlinear  problem,  we  offer  numeri¬ 
cal  evidence  that  there  is  asymptotic  improvement  in  the  rate  of  decrease  in 
the  error  as  a  function  of  computational  complexity.  For  linear  problems  we 
prove  that  a  version  of  our  method  converges  if  the  initial  data  is  sufficiently 
smooth. 

Our  method  is.  generally  speaking,  a  viscosity  method.  The  class  of  mono¬ 
tone  finite-difference  schemes  for  conservation  laws  are  also  viscosity  methods. 
Monotone  schemes  have  been  analyzed  by  Harten  et  al  [12],  Crandall  and  Majda 
[5],  Kuznetsov  [15],  Sanders  [19],  and  Lucier  [16].  These  schemes  converge  to 
the  entropy  weak  solution  of  the  conservation  law  (C),  as  formulated  by  Kruzkov 
[14].  Kuznetsov  provided  a  general  theory  of  approximation  for  approximate 

solutions  of  (C).  He  used  this  theory  to  provide  error  estimates  for  monotone 
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difference  methods  for  problems  in  arbitrary  space  dimensions  using  uniform 
grids  and  local  difference  operators.  His  techniques  were  used  by  Sanders  and 
by  Lucier  to  provide  error  estimates  for  difference  schemes  with  nonuniform 
meshes  and  nonlocal  difference  operators  respectively.  Sanders'  paper  provides 
a  lucid  treatment  of  Kuznetsov's  theory. 

One  of  the  considerations  in  developing  our  algorithm  was  that  it  must  exhi¬ 
bit  nonlinear  stability  properties  similar  to  those  of  monotone  finite-difference 
schemes.  Our  approach  may  be  conceptualized  as  follows.  We  take  a  uniform 
mesh  in  [a,6]x[0.f]  with  a  mesh  spacing  of  A f .  Meshpoints  are  removed  from 
the  mesh  where  they  are  not  needed  to  achieve  the  required  accuracy.  A  stan¬ 
dard  finite-difference  operator  is  used  to  advance  the  approximate  solution  from 
one  timestep  to  the  next.  (Special  techniques  are  used  when  the  meshes  differ 
from  one  timestep  to  the  next.)  Our  method  differs  from  previous  ones  in  the 
algorithm  for  choosing  the  mesh,  the  interpretation  of  the  approximate  solution, 
and  our  specific  choice  of  finite-difference  operator. 

Other  adaptive  methods  have  been  devised  for  evolution  equations.  Davis 
and  Flaherty's  algorithm  [6]  for  the  solution  of  evolution  equations  is  designed 
to  solve  smooth  problems  without  shocks.  They  use  an  L 8  analysis  to  choose 
their  mesh  and  to  provide  an  error  analysis  for  some  problems.  Our  algorithm  is 
similar  to  theirs,  in  that  it  tries  to  equidistrlbute  a  measure  of  the  error  among 
the  meshpoints.  However,  we  use  an  L 1  analysis,  since  the  solutions  of  (C)  are 
stable  in  Ll,  but  not  Lz. 

Other  general  algorithms  for  evolution  equations  were  proposed  by  Miiler 
[17]  and  Dupont  [9].  Dupont  supplies  convergence  analyses  for  his  methods. 
These  algorithms  are  mainly  finite  element  algorithms  that  converge  well  for 
solutions  that  ore  stable  in  L 8.  Since  the  solutions  of  (C)  are  stable  in  i,1  it  is  not 
immediately  clear  whether  such  algorithms  are  useful  for  approximating  (C). 
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Gannon  [11]  introduced  an  adaptive  finite-element  method  for  parabolic 
differential  equations  based  on  theory  for  elliptic  equations. 

Our  algorithm  was  also  motivated  by  the  work  of  Sanders  [19],  Douglas  [7], 
and  Douglas  and  Wheeler  [8]  on  monotone  finite-difference  schemes  with  nonuni¬ 
form  grids.  Sanders  and  Douglas  provide  convergence  results  for  such  methods 
with  a  fixed  grid,  while  Sanders  gives  an  error  estimate  of  0((A  +  At)*),  where  h 
is  the  largest  meshsize  and  A t  is  the  timestep.  Douglas  and  Wheeler  introduce 
an  algorithm  that  uses  grids  that  may  change  from  one  timestep  to  the  next,  a 
true  adaptive  mesh  method.  They  prove  that  the  solutions  of  their  algorithm 
converge  to  the  entropy  solution  of  the  conservation  law.  They  do  not  provide  an 
error  estimate.  We  compare  their  method  with  ours  in  the  final  section. 

When  Sanders,  Douglas,  and  Douglas  and  Wheeler  considered  a  nonuniform 
mesh,  they  interpreted  their  numerical  solutions  as  piecewise  constant  in  z,  and 
they  used  a  conservative  finite-difference  operator  that  is,  in  general,  incon¬ 
sistent  everywhere.  As  a  consequence,  maxfo -**_,)/ At  must  be  bounded  to 
achieve  stability  in  time  in  [8].  Since  our  method  can  have  arbitrarily  large  spa¬ 
tial  increments,  depending  on  the  smoothness  of  the  solution  and  the  nonlinear¬ 
ity  of  / ,  we  made  a  different  choice.  Our  method  interprets  the  solution  of  the 
numerical  problem  as  a  piecewise  linear  function,  to  take  advantage  of  at  least 
some  smoothness  in  the  solution.  We  also  use  a  consistent,  but  nonconservative, 
difference  scheme.  (Because  of  the  details  of  the  scheme,  it  is  still  conservative 
near  shocks.)  Section  4  contains  partial  results  that  bound  the  mass  error  in  a 
reasonable  way. 

Adaptive  numerical  methods  for  hyperbolic  conservation  laws  have  previ¬ 
ously  been  considered  by  Oliger  [10J  and  his  students,  iledatrom  and  Kodnguc 
[13]  is  a  survey  of  some  of  the  techniques  that  arc  used.  Dolstad  [2]  presents  a 
framework  for  methods  in  ouc  space  dimension.  His  schemes  incorporate 
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locally  varying,  recursively  defined,  space  and  time  increments.  He  uses 
Richardson  extrapolation  to  estimate  the  local  truncation  error  of  the  finite- 
difierence  scheme.  The  estimated  truncation  error  determines  the  local  grid 
size.  Berger  [1]  extends  Bolstad's  work  to  two  dimensions.  Among  other  things, 
she  deals  with  the  strictly  two-dimensional  problems  of  shock  capturing,  subgrid 
orientation,  and  overlapping  grids. 

The  stability  analyses  in  these  papers  are  L 2  analyses,  and  their  motivation 
seems  to  be  the  accurate  approximation  of  smooth  solutions  of  systems  of  linear 
conservation  laws.  Our  motivation  is  the  solution  of  nonlinear  problems  with 
shocks,  and,  in  dealing  first  with  the  scalar  equation,  we  use  a  Ll  analysis. 
Instead  of  finding  a  general  framework  for  the  such  methods,  we  use  a  specific 
finite-dilference  operator  and  mesh  selection  criteria. 

Oliger  and  his  students  employ  locally  varying  timesteps  as  well  as  spatial 
mesh  increments;  we  do  not.  We  could  have  used  locally  timesteps,  but  we  were 
not  able  to  prove  stability  and  convergence  results  for  these  methods.  Since 
asymptotic  improvement  in  convergence  rates  can  be  exhibited  with  fixed 
(small)  timesteps,  local  timesteps  were  not  used.  When  locally  varying 
timesteps  are  used,  our  algorithm's  implementation  is  very  close  to  Bolstad's. 

The  rest  of  the  paper  is  as  follows.  Notation  and  preliminary  results  com¬ 
plete  this  section.  Section  2  briefly  describes  the  finite-difference  operator  used 
here.  Section  3  presents  the  mesh  selection  algorithm,  and  proves  certain  use¬ 
ful  properties  about  the  resulting  mesh.  Section  4  contains  proofs  of  the  non¬ 
linear  stability  of  the  algorithm.  Section  5  shows  that  a  variant  of  our  methods 
converges  for  solutions  of  linear  problems.  Section  6  details  our  implementa¬ 
tion  of  the  algorithm,  and  Section  7  describes  our  computational  results. 

The  BV  seminorm  of  a  function  u  is  defined  as  |u  =  f  |u'(x)  |  dx. 

a 

where  the  integrand  is  interpreted  as  a  finite  measure.  The  set  of  all  such 
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functions  Is  denoted  by  BV(R).  If  it  is  in  BV{ R),  then  there  exist  two  bounded 
functions,  u *  and  u~  such  that  u  =u*  +  u~  and  u*  is  nondecreasing.  u~  is  nonin¬ 
creasing.  We  define  u *  =u+  -u~,  the  total  variation  function  of  u . 

Throughout  the  paper  we  assume  the  normalization  that  ||/*|li«»  — 1>  This 
can  always  be  achieved  with  a  change  of  the  time  scale,  and  is  used  only  for  con¬ 
venience  in  stating  stability  conditions. 

2.  The  Unite-difference  Scheme. 

We  use  a  standard  upwind-difference  scheme  to  advance  the  approximate 
solution  from  time  tn  to  fn+1.  We  are  given  a  suitable  mesh,  chosen  by  the  rules 
in  the  next  section,  to  represent  the  solution  at  time  tn .  It  is  characterized  by 
the  meshpoints  x{*  and  the  values  of  the  approximate  solution  U?  at  those  mesh- 
points.  We  interpret  these  points  as  determining  a  continuous  piecewise  linear 
approximation  to  u(x,fn).  A  estimate  of  the  solution  at  time  fn+1  is  calculated 
by 

ur'-v? .  /*w)-/'W-,>  . 

~a  v  wl  (  ’ 

* 

for  all  x f  (except  the  two  endpoints  of  the  interval).  We  have  decomposed  /  into 
its  increasing  (/+)  and  decreasing  (/")  parts.  If  /  is  monotone  increasing  or 
decreasing  then  this  method  is  an  upwind  difference  method.  A  similar  scheme 
has  been  used  by  Engquist  and  0sher[10j. 

The  linear  intcrpolant  of  the  values  t/"+1  at  the  points  x"  is  a  function  we 
call  uH.  The  mesh  selection  algorithm  of  the  next  section,  when  applied  to  uh, 
givos  us  a  new  mesh  xtn+1  and  function  .''dues  U?*1  for  the  approximate  solution 
of  (C)  at  time  fn+l. 

This  process  is  repeated  until  fn+1  *  T. 
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3.  The  Heah  Selection  Algorithm. 

This  section  describes  our  mesh  selection  algorithm.  In  our  implementa¬ 
tion,  the  mesh  at  time  f**1  is  built  from  the  mesh  at  time  tn,  but  the  method  of 
approximation  is  general  and  applies  to  any  function  with  reference  to  a  time¬ 
stepping  procedure.  We  present  the  algorithm  here  in  its  general  form. 

Our  method  of  mesh  selection  is  similar  to  well  known  ones  for  adaptive 
linear  approximation^].  The  mesh  approximately  equidistributes  an  estimate  of 
part  of  the  error  incurred  by  the  finite  difference  scheme,  thus  following 
methods  used  in  static  problems  [4]  and  other  evolution  equations[6].  The 
method  presented  here  was  designed  for  problems  with  shocks,  while  previous 
methods  were  designed  for  possibly  nonlinear  problems  with  smooth  solutions. 

The  problems  in  [6]  succumbed  to  an  L 2  analysis,  while  the  conservation 
laws  that  are  the  target  of  this  method  are  stable  only  in  L1.  The  mesh  selec¬ 
tion  algorithm  therefore  chooses  a  mesh  that  is  "right”  for  Ll.  Our  specific 
choice  of  the  mesh  will  allow  us  to  prove  the  stability  results  of  section  4.  an 
important  goal. 

Let  u  be  any  bounded  function  defined  on  [a  ,6aj  that  is  constant  outside  of 
[a,b],  and  let  e  be  a  small  parameter.  Choose  the  mesh  according  to  the  follow¬ 
ing  algorithm: 

ALGOEHTHU  H:  This  algorithm  chooses  meshpoints  at  which  to  approximate  an 
arbitrary  function u  defined  on  [a,6], 

1.  The  meshpoints  consist  only  of  the  points  a  and  b  and  the  centers  of  admis¬ 
sible  intervals.  Admissible  intervals  are  defined  by  (2)  and  (3)  below. 

2.  The  interval  [o ,6]  is  an  admissible  interval. 

3.  For  any  admissible  interval  /,  let  3/={z  | dist(x,/)=ml'|;c-v/ 1 <  |y  j  j.  If  1  / 

l/c/ 


and 
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|/|/l“*»l  +  (3.1) 

a i 

then  the  left  and  right  halves  of  I  are  admissible  intervals.  The  above 
integral  is  finite  if  u„  is  a  finite  measure;  it  is  to  be  interpreted  as  <»  other¬ 
wise.  Note  that  3/  is  an  open  interval. 

When  this  algorithm  is  used  for  our  adaptive  method,  Af=e/  4,  so  that 
At  <min(xj* -*T-i  )• 

A  minimal  interval  is  an  admissible  interval  that  contains  no  proper  admis¬ 
sible  subintervals.  It  is  clear  that  the  width  of  any  given  admissible  interval  is 
2** (6 -a)  for  some  nonnegative  k.  Also,  ]/|  >  c/2.  It  follows  that  the  mesh  is  a 
subset  of 

Sc  =  Ja+m2~fc(b-a)  |  m  =0, 1,  •  •  •  ,2*  1and(6-a)2_*ae/4>(6-a)2"k''1j 

Some  lemmas  about  the  structure  of  the  mesh  generated  by  this  algorithm  fol¬ 
low. 


Leuua  3.1.  If  A  and  B  are  two  adjacent  minimal  intervals,  then 


L*cl£L 

2  \B\ 


<2. 


Proof.  We  let  \B  |  >2  |A|,  and  derive  a  contradiction.  Since  the  width  of 
any  admissible  interval  is  (b-a)2~*  for  some  0,  |1?|>4|A|.  Consider  now 
the  admissible  interval  from  which  A  was  derived  by  Step  3,  and  call  it  C.  Since 
C  was  divided  into  two  admissible  subintervals,  the  integral  in  Step  3  is  greater 
than  e.  However,  \B  |  a 2  |  C|  and  3Cc327,  so  that  the  corresponding  integral  for 
B  must  also  be  greater  that  e,  a  contradiction  to  the  minimality  of  B. 


/// 


Except  for  the  two  points  a  and  b,  the  set  of  meshpoints  has  the  natural 
structure  of  a  trer.  The  point  fa+b)/  2  is  the  root  of  the  true.  You  con  also  think 
of  the  interval  [a.pi  as  ie  root  of  the  tree.  (Since  there  is  a  one-to-one 
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correspondence  between  the  meshpoints  and  the  set  jf  admissible  intervals,  we 
will  describe  the  structure  of  the  mesh  equivalently  in  terms  of  intervals  or 
meshpoints.)  If  an  interval  /  is  divided  into  two  admissible  subintervals  by  Step 
3,  then  these  two  intervals  are  the  left  and  right  children  of  /,  and  /  is  their 
parent.  A  meshpoint  with  no  children  (the  center  of  a  minimal  interval)  is  a 
leaf. 

Lemma  3.2.  Every  second  meshpoint  is  a  leaf. 

Proof.  The  statement  of  the  lemma  is  true  after  Step  2,  and  it  is  left  invari¬ 
ant  by  Step  3. 

/// 

This  implies  the  obvious  result  that  there  is  a  unique  covering  of  [o.6]  with 
minimal  intervals. 

Lemma  3.3.  If  the  set  of  meshpoints  is  {*<(  with  A*  =  x*  then 

l/2^Ai/Ai+1<2. 

Proof.  If  xt  is  a  leaf,  then  A*  =hi+l.  If  x^  and  xi+1  are  leaves,  then  the 
result  follows  from  Lemma  3.1. 

/// 

The  linear  interpolant  uf  of  u  is  defined  by  requiring  that  ue(x*)  =  ti(x<)  for 
all  Xi,  and  that  ue  be  a  linear  function  between  the  meshpoints  so  that  ut  is  con¬ 
tinuous  on  [a  ,6].  The  function  u(  has  the  following  approximation  properties 
(cf.  [3J). 

Lemma  3.4. 

Ilu -’ucIIl1([o.&])  ~  ~a  I  +  S  ilu  -u*lli‘(/)  (3-2) 

I  minimal 
|/|<« 

Proof.  If  /  =  (xt_,,x1+I)  is  minimal  and  |/|£c  then  Step  3  implies  that 
|/|/l«*»  \dx<,t  The  LX{I)  error  for  linear  interpolation  at  the  points  x<_|.  xi, 

i 


and  xt*!  is  bounded  by 


Summing  over  all  I  gives  the  first  term  in  the  expression.  The  second  term  con¬ 
tains  all  the  intervals  not  yet  considered. 

/// 


There  are  two  interesting  cases  when  the  second  term  is  known  to  be  small. 

Lemma  3.5. 

(a)  If  u  is  a  continuous,  piecewise  linear  function  such  that  ux  is  discon¬ 
tinuous  only  at  the  points  in  Sc,  then  the  second  term  in  (3.2)  is  zero. 

(b)  If  u  is  in  BV(R)  then  the  second  term,  in  (3.2)  can  be  bounded  by 

1  /  2 1 |it  1 1/? v(i^ e ■  (Here  we  assume,  without  loss  of  generality,  that 
z*-h 

u.(x)  =  lim— -  f  u(t)dt  for  allx.) 
a-*o  <th 

Proof. 

(a)  Since  u  is  linear  on  each  half  of  the  minimal  intervals  /  with  |/|  <e. 
ut  =  u  on  these  intervals. 

(b)  It  is  an  exercise  to  show  that  \\ut-u\\Lx^<t\/ 2t\\u\\Bv^jy  Since  u  is  in 
BV( R),  we  may  add  these  individual  bounds  to  obtain  the  lemma. 

/// 


A  few  other  stability  properties  will  be  used  in  the  sequel. 

Lemma  3.6.  ||Welbv([«.6l)-llwllz?va«.6]).  and  ||ue||i.(lo  b])<  6]). 


Proof.  It  is  clear  that  linear  interpolation  has  these  properties. 


f 
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4.  Stability  Properties. 

The  numerical  method  presented  here  has  several  stability  properties  that 
mimic  those  for  conservation  lavrs.  In  brief,  solutions  of  the  numerical  method 
satisfy  a  maximum  principle,  are  total  variation  diminishing,  and  are  stable  in 
time.  Although  the  numerical  scheme  is  not  conservative,  one  can  bound  the 
mass  error  for  most  problems. 

Boundary  effects  will  be  analyzed  in  the  following  way.  Outside  of  [a, 6]  the 
mesh  will  be  extended  so  that  all  the  mesh  intervals  to  the  left  of  a  sure  of  the 
same  width  as  the  mesh  interval  immediately  to  the  right  of  a,  denoted  by  ha-  A 
similar  extension  will  be  made  to  the  right  of  b  .  We  assume  that  the  function  l/tn 
is  constant  to  the  right  and  to  the  left  of  [a, 6].  The  mapping  from  V f  to  l/"*1  is 
now  divided  into  three  parts:  the  finite-difference  scheme  (2.1)  transforms  1/"  on 
the  extended  mesh  to  U?*l\  the  values  at  the  points  a  and  6  are  reset  to  their 
original  values;  and  the  remeshing  procedure,  applied  to  the  mesh  values  in 
[a. 6  J,  yields  C/"*1. 

For  the  method  to  work  properly,  some  criteria  is  needed  to  decide  when  a 
shock  or  other  disturbance  is  getting  too  close  to  the  boundary  of  the  interval. 
Throughout,  we  will  assume  that  the  minimal  intervals  adjacent  to  the  boundary 
points  a  and  b  have  diameters  greater  than  e.  This  is  a  simple  and  effective  cri¬ 
teria.  If  this  criteria  is  in  danger  of  being  violated,  the  interval  [a.b]  is  to  be 
enlarged. 

We  will  follow  the  development  of  [0]  for  many  of  our  pr  oofs. 

The  following  simple  lemma  will  have  important  applications. 

LEMMA  4. 1 .  If  g  is  f  f~  or  / 1 ,  and  u  is  the  linear  irdcrpalani  of  the  points 


(x«n, Vi)  then 
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gUgliWW)  g(u?)-g(UT-i) 

Afci  "  fkn 


*JVi 

/  I«W  I  +  !/"(«)!  ufdx  (4.1) 
*f-i 


Note  that  this  is  multiple  of  the  quantity  that  we  use  as  a  subdivision  criteria. 
Proof.  Since  g  has  a  bounded,  piecewise  continuous  second  derivative, 


( IP  —  IP'S 

g{U?*i)  =  gW)  +  fti"  1  g  W)  -  '  —4  /  (**"♦, -x)p(u(x))«  dr. 

Now,  g'(u)  is  bounded  by  1,  and  for  xe(xiB,x<V i),  g  (u)„  =  g"(u)u^,  since  =  0. 
Also,  |g*’|  is  either  |/"  |  or  0.  The  previous  equation  and  a  similar  one  for 
g  ( U^-i )  can  be  rearranged  and  summed  to  yield 


g(v?+>)-9(u?)  g(ur)-9{ur- 1) 

AJVi  '  V 


/  |U»I  +  l/'Wl^dx. 


Here  we  have  expressed  the  difference  of  the  left  and  right  derivatives  of  u  at  x? 
as  the  integral  of  the  second  derivative  “delta”  measure. 


/// 


The  sharper  bound 


/W«)-/W) 

/WWW- 1) 

4. 

/w*iww) 

/WWW-i) 

T 

v+i 

<  f  \u„\  +  |/"(u)|ul3dx  .  (4. 2) 

may  be  derived  by  noting  that  /  =  f*+f~>  and  if  the  first  or  second  derivative  of 
f*(u(x))  is  nonzero  at  x,  then  the  first  and  second  derivatives  of  /~(u(x))  are 
zero  there. 

Theorem  4.1.  For  alln>Q,  sup £/”£  sup U°. 


Proof.  By  Equation  (2.1), 
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Since  f*  is  increasing,  /  "is  decreasing,  and  At/h?£  1,  the  preced¬ 

ing  equation  yields 

. ...  >.\f+(u?)-f*wr-ivu?vu?„)  .  riur-^urvur^-f-iu?) 

U?  <U?-At  - — - + - — 

a*  'k  +  i 

<  U?  +  /*(  U?-i  V  £//V  L/JV , )  - / ‘  ( i/i*) 

<  Ui-i  V  C/?*V  LWi- 


Therefore,  the  finite-difference  scheme  satisfies  a  maximum  principle.  Reset¬ 
ting  the  values  of  t/"  at  a  and  b  does  not  violate  a  maximum  principle.  Lemma 
3.6  states  that  the  subsequent  transformation  to  IP*1  satisfies  a  maximum  prin¬ 
ciple. 


/// 


Theorem  4.2.  ForaUn>Q,  ||C/nlbv«)£l|i/0llRv(«- 
Proof.  If  we  let  S/>  =  / *(  U?)  -/  ♦( U?- , ).  then 


-  UT-V  =  v?-  VlLi -At 


W 

8/t-i 

Sfili 

a/f 

h? 

A?-i 

th\i 

\n 

Because  f*  and  /  are  monotone,  ft=f*  -/",  ||/,||2,5>  —  1  and  At //i/*<  1,  we  can 
take  absolute  values,  and  sum: 


Zwr^-vr-Yi  £ZWin-vr-i  -^[/‘W)-/W-i)]|+a*e  i 

i  i  i 

=  E\  vr-uz.il. 


tifi*  1 
fcJVi 


t 

Thus  the  mapping  U*-*  IP  is  total  variation  diminishing.  Resetting  the  values  at 
the  endpoints  does  not  increase  the  total  variation.  Since  the  remeshing  pro¬ 
cess  is  variation  diminishing,  by  Lemma  3.6,  the  theorem  is  proved. 


/// 


THEOREM  4.3.  || Un*1  -  t/n||ii{[Bit,J)^(3/2||(;0||i7K(B)4(b  -a)/ 0)A t  +  2A tt. 


Proof.  First,  since  the  mesh  is  graded. 

l!£/"7r-  tf"lli>(t..n)  *  I  ^nT"  W\ 

£  4f-£  |/+(t^)-/+(WT-i)l(i+^+l/'(«?Vi)-/-(«r)l(i 

^  i  "i 

*  ^Li;i/‘(t/in)-/‘(^-i)i 

When  the  boundary  values  are  restored,  we  commit  an  error  at  the  left  end¬ 
point  of  at  most  At  |  U?  -  l/fl  |.  Under  the  assumption  that  the  width  of  the 
minimal  interval  adjacent  to  a  is  bigger  than  e,  Step  3  of  the  mesh  construction 
allows  us  to  bound  the  error  by  Ate. 

Finally,  we  may  apply  Lemmas  3.4  and  3.5  to  yield 

II  £/"+1  -  (6  -  a)  At  /  6. 

/// 

Because  a  consistent  but  nonconservative  scheme  is  used  for  the  time  stepping, 
there  will,  in  general,  be  some  mass  balance  error.  This  error  can  be  bounded 
by  the  following  theorem. 

THEOREM  4.4.  .Assume  [a, 6  ]  =  \jSj  (J  \jR}  where  m- 1  disjoint  regions  R}  are 

i  i 

made  up  of  minimal  intervals  of  width  less  that  c,  and  are  surrounded  by  m 

regions  Sj  that  are  covered  by  minimal  intervals  whose  width  is  greater  than  e. 

Assume  that  there  are  no  more  than  k  minimal  intervals  in  [jSj.  Then 

; 

j  f  ^rn-_  l/«  dx  -  At  (/  (u  (a ))  -/  (u  (b )))  I  25  At  (2  k  t  +  (2  m  e  ||  l/n  ||flK(B))»)  (4.3) 

B 

Proof.  Using  the  notation  of  Theorem  4.2,  we  have 
jrfU""(zj-U*(x)dx 


.  14- 

^  V  A&!  n  2  }’ 

=  k/*(u(«))  +  /-(«(a))-Cr(u(6))  +  /-(ii(b)))]-S(^*Ri  +  ^LV) 

*  t  "t  rtt  +  l 

=  kr<  «(«))-/<«(6)))-e^Wm-E^^v. 

«  i  A*  i  rtt  +  1 

We  will  show  that  the  first  sum  minus  l/2(/+(u(a))-/*(u(6)))  can  be  bounded 
as  in  the  statement  of  the  theorem. 


Pick  a  particular  Sj  with  left  and  right  endpoints  xt  and  Xy  and  consider  the 

sum 


*T 

£ 

*r=*« 


8/i* 


NVi- 


(4.4) 


Since  every  minimal  interval  contained  in  St  is  more  than  e  wide,  Step  3  of  Algo¬ 
rithm  M  and  Lemma  3.5  show  that  (4.4)  is  equal  to 


ft  V&i 

*f=*, 


iu  v 


*t+jw  . 


where  |  A/  |  is  less  than  e  (r  -i  + 1). 


Any  interval  Rj  consists  of  minimal  intervals  of  width  less  than  e.  Perforce, 
all  the  intervals  must  have  the  same  width.  Thus,  if  zt  and  are  the  left  and 
right  boundaries  of  R{, 


*r— 1 


E 


*/*=*<♦! 


*r-! 

E 


* i*i 


Matching  these  results  gives 


£ 


E 


+  A/  + 


si 


E  (Vr*  l 
=  (*!•«,) 


-3/.+). 


(4.5) 


Here,  \M\^2kt,  since  I JS}  is  covered  by  k  minimal  intervals.  The  remark  fol¬ 
lowing  Lemma  4.1  implies  that  the  same  bound  holds  when  the  terms  incorporat¬ 
ing  J~  are  also  included.  The  first  sum  on  the  right  of  (4.5)  is  equal  to 
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a)). 

if*  6/.* 

For  the  interval  S*.  assume  |  — - — — |  =  Kjt  for  some  positive  K>. 

*  V 

Then,  without  loss  of  generality,  we  may  assume  that  Lemma  4.1 

"T  +  l  * 

iff.  iff 

and  Step  3  of  the  mesh  construction  imply  that  |  — - — ^ 5 1  if  x*£Sj. 


ui 


Also,  e /  2  for  xfe-Sj .  Therefore, 


Jtj/8 


i-0 


2s  e  Kf/  8 . 


Thus,  the  last  sum  in  (4.5)  is  bounded  by 

e/2 £  (2  cm)M(e/ &£*/)» 
i  i 

S{Ztm.\\f*{U»)\\Bvmy 

Since  ||/+(C/t*)||bk{b)  +  !l/-(£/n)llsy(B):S  ||C/"||*vnp.  the  theorem  is  proved. 

/// 

It  can  be  shown  that,  for  piecewise  smooth  functions,  Algorithm  M  generates 
fewer  than  Kz~%  meshpoints,  for  some  K  (cf.  de  Boor[3j).  If  this  bound  is 
assumed  throughout  the  calculation  of  u(,  Theorem  4.4  bounds  the  mass  error 
at  a  fixed  time  T  by  a  multiple  of  e*  or  it*.  The  expected  rate  of  convergence  of 
the  method  is  0(Af#). 

Unfortunately,  we  have  not  been  able  to  calculate  an  a  priori  bound  on  the 
number  of  meshpoints  in  a  mesh  generated  by  Algorithm  M.  If  the  function  ue  is 
rough  enough,  c-1  meshpoints  may  be  needed  to  approximate  it.  The  numbers  k 
and  m  can  be  calculated  during  the  course  of  the  algorithm  in  a  negligible 
amount  of  computer  time,  however,  and  (4.3)  may  be  used  as  an  a  posteriori 


bound. 
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Theorem  4.4  does  not  consider  the  mass  error  generated  by  the  adjustment 
of  the  function  values  at  the  points  a  and  6 .  Theorem  4.3  shows  that  this  error 
is  O(At)  if  the  minimal  intervals  next  to  a  and  6  have  width  greater  than  e. 

Estimate  (4.3)  does  not  include  the  mass  error  caused  by  the  remeshing 
process.  Such  an  error  occurs  when  admissible  intervals  in  the  mesh  at  time  tn 
are  no  longer  needed,  and  do  not  appear  in  the  mesh  at  time  tn+1.  Lemma  3.4  is 
sharp  for  any  given  time  step;  the  mass  error  for  a  time  step  of  size  At  can  be 
comparable  to  At.  Using  this  bound  over  the  interval  [0,7']  gives  an  estimate  of 
an  0(1)  mass  error,  a  totally  unacceptable  result. 

Computational  experience  with  piecewise  smooth  solutions  of  the 
differential  equation  has  shown  that  over  the  interval  [0,7*]  there  is  a  constant  C 
such  that  fewer  than  C2*  intervals  of  width  2"*(6-a)  are  removed  from  the 
mesh  in  [0.7’].  That  is,  a  minimal  interval  is  not  subsumed  into  a  larger  minimal 
interval  for  a  period  of  time  proportial  to  the  its  width.  If  this  property  holds,  an 
argument  similar  to  that  of  Lemma  3.4  shows  that  the  mass  error  due  to  mesh 
changes  on  [0,  T]  is  bounded  by  CAflog(Af“l)  for  some  C.  We  therefore  have  the 
following  theorem. 

THEOREM  4.5.  For  an  interval  [O,?*],  assume  that  there  are  constants  C,.  C2. 
and  Ca  such  that  for  any  e>0 

(a)  there  are  no  more  than  minimal  intervals  of  width  bigger  than  e  for 

any  tn  in  [0,7’], 

(b)  there  are  never  more  than  Cz  disjoint  regions  covered  by  minimal  intervals 
of  width  less  than  e,  and 

(c)  at  most  C3 2*  admissible  intervals  of  width  (6  -a)2~*  are  removed  from  the 
mesh  in  [0, 7’]. 


'then  there  exists  a  constant  C,  depending  on  C\,  C3,  C3.  and  ||u such  that 
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tha  mass  srror  of  the  numerical  approximation  (2.1)  can  be  bounded  by  Ce*. 

/// 

The  conditions  (a),  (b),  and  (c)  can  easily  be  checked  a  posteriori. 

5.  Convergence  for  linear  Problems. 

In  this  section  we  prove  that  solutions  of  the  linear  problem 

Uf  +  cru,  =  0,  xeROG,  (5.1) 

u(x,0)  =  u0(x),  zeR 

are  approximated  well  by  variants  of  our  adaptive  mesh  algorithm.  If  the  initial 
data  is  in  BV(R),  a  posteriori  bounds  will  be  given  on  the  error,  and  if  the  initial 

dUn 

data  is  slightly  smoother,  with  — — e.BV{ R),  a  priori  bounds  are  available.  For- 

OX 

mol  justification  of  the  mesh  selection  criteria  used  in  Algorithm  M  is  also  given. 
Consider  first  when  uaeBV(R).  Our  algorithm  is  as  follows: 

1.  Pick  c>0,  and  let  u$  =  "here  yt(x)  is  t~l/z  when  |x  |Sel/z/2  and  0 

otherwise. 

2.  Using  uj  as  the  initial  data,  follow  the  algorithm  in  Section  2  to  find  «*. 
This  will  be  our  approximation  to  u . 

It  is  easily  shown  that  ll^^l*K»*e~1/8ll“oll*v».  and 
Ulig  —  Ujf 

The  following  theorem  applies  to  such  problems. 

Theorem  5.1.  If,  at  each  time  step,  there  are  never  more  than  C1e~l/Z  minimal 
■intervals  of  width  greater  than  e;  and  if  the  total  ^‘([a.b])  error  due  to  mesh 
changes  is  0(eui),  then 

IMO  +  l)el/z(l|ttolleviB) +  1)  (5.2) 

thirthennore ,  the  computational  complexity  of  the  schema  is  0(e_5/8) 
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Theorem  5.1  gives  a  posteriori  estimates  of  the  error  and  computational 
complexity  of  the  scheme.  Our  computational  experience  suggests  that  the  con¬ 
ditions  of  the  theorem  are  always  satisfied.  We  suspect  that  the  structure  of 
solutions  of  the  conservation  laws  ensures  this.  The  computational  bound  is 
better  than  for  standard  monotone  methods;  for  these  methods  the  error  is 
0(Af l/z)  but  the  complexity  is  0(Af  _z). 

The  proof  of  this  theorem  relies  on  a  number  of  lemmas. 

Leioia  5.1.  //  u,ei?K(R)  and  v  is  the  linear  interpolant  of  u  on  some  mesh, 
{xj.  then  |(u,|bv®S||u,||^|W. 

1  * 

Proof.  Since  vx(x)= - f  ua(t)dt  for  xe(x1_llxt)1  the  result  follows 

x*  “*<-»  *_t 

immediately. 

/// 


Leioia  5.2.  If  U11  is  generated  by  Algorithm  U  and  (2.1),  then 

II^xILbk(b)  ^  ll£4°llflKw- 


Proof.  Since  (5.1)  is  linear,  VJ*  = 


_  ur-ur-x 


satisfies  the  same  difference 


equation  as  V f.  Hence  ||V*1|Ij?kws  llVMIavro-  By  the  previous  lemma, 

The  result  then  follows  by  induction. 


/// 


Our  analysis  views  upwind  finite-difference  schemes  for  conservation  laws  in 
a  seemingly  new  way.  In  particular,  if  U?  and  are  interpreted  as  piecewise 
linear  functions,  then  the  equation  (2.1)  defining  U?+l  is  equivalent  to  the  follow¬ 
ing  algorithm: 

To  obtain  U^*1,  shift  U?  to  the  right  by  aA t  (if  a  is  positive)  and  interpolate 
the  shifted  function  at  the  points  x(n. 
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Note  that  error  occurs  only  In  the  interpolation  process;  there  is  no  error  in  the 
values  of  the  U?*1  themselves  if  the  values  of  the  U?  are  exact.  This  calculation 
is  illustrated  in  Figure  1. 


Proof  of  Thaorem  5.1.  We  first  define  a(w)(x)  =  v(*-oAf).  If 

e»  =  ||u(«AO  ~  then 

=  ||a(u((n-l)A«))-a(cr»-1)||4,(lMn  +  ||a(l^-*)  -  ^llxi([..b]) 

=  en-l  +  3n 

where  5n  is  the  local  error  110(1/"  1)  —  U*  llxi([»,>])-  define 

5Ui-{Ux~  Ui-\)/hf.  Now,  from  Figure  1, 

„t^-o<t/n)!liI([...I)  =  EaA«W(l““-)l«Kr-«wr-t  I  <5  3) 

*1*  ^ 

S2aA t*  £ 

Hf<X6t 

+  «Af  £  Ai‘W-W~il 

KptZM 

By  using  Lemmas  5. 1  and  5.2,  one  can  bound  the  first  term  by 

CA*e||t/"||tf»TO  £  Ct  1/*Hu0!tay(s»A* 
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Because  of  Step  3  of  Algorithm  M,  and  our  assumption  about  the  number  of 
intervals  of  width  greater  than  e,  the  second  term  can  be  bounded  by 

aAf£xC£-1/2  =  Ce1/2Af 

Thus  ||t^-a(^)||£l{[B  >])^CAt£1/a(j|u0|bnB)+  1). 

Summing  this  expression  over  n<t/ht  gives  one  term  of  (5.2).  We  have 
assumed  that  £||t/n+1  -  ^n+1llii([ai6])- Ce54.  As  stated  previously. 

n 

||u0-uS||<£1/2|fuolk?V(B)'  Lemmas  3.4  and  3.5  show  that 
\\U° -ug \\Li^a  b])£ t(b  -  a)/ 32+  l/2£||u0|bv(H).  The  error  bound  is  proved. 

The  complexity  of  the  scheme  is  0(&t~lN)  where  N  is  the  maximum 

number  of  meshpoints  in  any  mesh  {a:"),  since  the  amount  of  work  is  linear  in 

the  number  of  meshpoints.  By  hypothesis,  there  are  fewer  than  Ce~1/2  minimal 

intervals  of  width  greater  than  c.  Furthermore,  since  f  \u„  |  dx>l  if  /  <e,  there 

3/ 

must  be  fewer  than  3HU^llByW£3c~i/2llu0llBY^  minimal  intervals  of  width  less 
than  c.  Because  every  second  meshpoint  is  the  center  of  a  minimal  interval,  the 
theorem  is  proved. 

/// 

We  now  consider  the  case  when  u  is  smoother.  Let 

Ll  2(R)  =  jueZ,‘(R)rW(R)  I  u,e£F(R)}. 

We  define  here  a  modification  of  Algorithm  M  for  functions  in  Z.l,z(R). 

Algorithm  M' 

1.  The  points  o  and  6  are  meshpoints.  The  center  of  every  admissible  interval 
(to  be  defined  below)  is  a  meshpoint. 

2.  The  interval  [a.b]  is  an  admissible  interval. 

3.  If  /  is  an  admissible  interval,  \I  |S4Af ,  3/  =  \x  |  dist(x,/)< \I  |  j,  and 
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then  the  left  and  right  halves  of  /  are  admissible  intervals.  As  before,  u„  is 
to  be  interpreted  as  a  measure. 

This  new  mesh  algorithm  makes  the  smallest  admissible  intervals  the  same 
size  as  the  “average"  interval,  i.e.  when  is  bounded.  Since  ue£1,!l(R),  the 
very  small  spatial  and  temporal  mesh  increments  of  the  previous  method  are 
not  needed.  The  following  adaptive  algorithm  is  therefore  a  true  "smooth  solu¬ 
tion"  algorithm. 

Our  new  algorithm  is  as  follows: 

(a)  Given  u0(x).  And  xf  ising  Algorithm  M'  and  choose  Ut°  using  piecewise 
linear  interpolation  of  u0. 

(b)  For  each  timestep:  Given  meshpoints  {*tn|  and  function  values  { U?\  defined 
at  those  meshpoints,  calculate  U"+l  using  (2.1).  (Equation  (2.1)  reduces  to 
the  upwind  differencing  scheme  in  this  case.)  Interpreting  U?  as  a  continu¬ 
ous,  piecewise  linear  function,  use  algorithm  M'  to  find  a  new  mesh  {x"+1j 
for  Ui*1,  and  define  U"*1  on  that  mesh  by  linear  interpolation  of 

The  following  result  holds. 

Theorem  5.2.  Let  u0e£l1,8(R)  and  u(x,f )  be  the  solution  of  (5.1).  If  nAf  =  T, 
and  £/"  is  the  solution  of  the  adaptive  mesh  algorithm  above,  then 

IMO  -  *  3<r  +  DA* 

Proof  As  in  Theorem  5.1, 

!l^riT-a(t/n)||il([.ib]^2aAfa  £  1 5 U? - 5 £/<L,  |  +  aAf  £  h?  | 

The  first  term  is  bounded  by  2aAfa||«o'lbv(B)  because  of  Lemmas  5.1  and  5.2.  By 
using  Step  5  of  Algorithm  M',  the  second  term  can  be  bounded  by 


■-**■»*  Jp,‘l 


<*At[  £  Af|l/*f  E 


1/8 


S(U.<6-„y'««#%® 

S  oAi8||lio'lbK(B)' 


(6  -a)1'8 


Thus  ||£/n+l  -a(^)||Il([a<6])^ 3aAf8j|«o’l!i»^». 

Finally,  let  5  be  the  set  of  all  minimal  intervals  in  the  mesh  at  time  tn+1 
that  are  not  minimal  at  time  t".  Then,  because  of  Lemma  3.4  and  step  3  of  Algo¬ 
rithm  M'.  the  expression  for  the  error  in  the  trapezoid  rule  yields 
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Thus.  5n<3Atz||iio'[|flV(E-  After  calculating  the  error  caused  by  the  approxima¬ 
tion  of  u0,  the  theorem  follows  by  induction. 

/// 

The  same  analysis  can  be  used  formally  for  nonlinear  problems:  estimate 
the  difference,  now  nonzero,  between  the  values  C/tn+1  and  (t/")(xtn),  (Sw 
advances  the  solution  of  (C)  by  time  At )  and  then  separately  estimate  the  inter¬ 
polation  error  as  in  Figure  1.  The  error  in  C/tn+l  is  bounded  formally  by 

*T*t 

At  /  \f(Un(x,ln))xx  |  dx  +  0{hlz).  (0.4) 

*T-i 

The  L1  error  on  (ztn-i  ,x?M )  is  therefore  bounded  by  (hin  +  /itIVi)/2  times  (5.4), 
bounded  in  turn  by  CAt  times  the  integral  in  (3.1). 

Thus,  if  ||t/*  *  lllz/K(/)  —  ^W^xWovi^i)  h>r  ouch  miuuual  interval  /,  Llio  error- 
caused  by  the  nonlinearity  is  of  the  same  order  as  the  error  caused  by  the 
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dissipation  in  the  difference  scheme.  Thus,  our  mesh  selection  criteria  achieves 
a  balance  between  errors  caused  by  nonlinearity  and  dissipation. 

6.  Implementation  Details. 

Our  scheme  wets  implemented  using  the  programming  language  Pascal.  The 
resulting  programs  were  run  on  a  VAX  1 1  /780  with  a  floating  point  accelerator 
under  the  VMS  2.5  operating  system  and  Pascal  1.2  compiler,  and  under  the 
Berkeley  Unix  4.1  operating  system  with  its  Pascal  compiler.  In  this  section,  we 
describe  the  algorithms  and  data  structures  used  in  our  implementation.  We 
show  that  the  integrals  in  (3.1)  can  be  evaluated  exactly,  and  hence  that  the  sta¬ 
bility  results  of  the  previous  section  hold  without  a  separate  theory  describing 
the  effects  of  numerical  quadrature.  We  also  estimate  the  computational  com¬ 
plexity  of  the  scheme. 

The  following  steps  advance  the  approximate  solution  from  one  timestep  to 
the  next.  First,  the  finite-difference  formula  advances  uh  from  tn  to  fn  +  1  on  the 
mesh  |z/*j.  Secondly,  the  integrals  in  (3.1)  are  calculated  for  the  mesh  |zt”|.  In 
our  implementation,  the  integrals  are  first  calculated  over  the  (open)  interval  I 
instead  of  31,  and  the  integral  over  3/  is  constructed  when  it  is  needed.  Finally, 
the  mesh  {z<*+,{  is  constructed.  To  do  this,  the  union  of  the  meshes  jz^J  and 
Iz/1*1}  is  built  up  by  adding  the  appropriate  meshpoints  to  \x?\.  The  values  of 
the  integrals  (3.1)  are  derived  for  the  new  meshpoints  as  they  are  introduced.  A 
second  pass  is  made  to  remove  points  in  jztn(  that  are  not  needed  in  {z"+1{. 

Except  for  the  two  meshpoints  a  and  b ,  the  mesh  is  naturally  organized  as 
a  binary  tree.  This  is  because  it  is  defined  recursively  by  subdividing  admissible 
intervals  into  two  subintcrvals  at  each  step.  Each  admissible  interval  (or, 
equivalently,  the  meshpoint  at  its  center)  is  a  node  in  the  tree.  The  interval 
[a, 6]  is  the  root  of  the  tree.  If  an  admissible  interval  is  subdivided  into  two 
admissible  subintervals,  then  these  subintervals  arc  the  left  and  right  children 
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of  the  interval.  A  node  without  children  is  a  leaf .  and  corresponds  to  a  minimal 
interval.  Nodes  that  are  not  leaves  are  interior  nodes.  The  parent  of  an  interval 
is  the  admissible  interval  from  which  it  was  derived. 

We  first  describe  the  information  stored  in  a  node  corresponding  to  a  mesh- 
point  x,n  at  the  center  of  an  interval  /  =  (xt,xr}.  This  information  consists  mainly 
of  pointer  variables  and  numeric  variables.  A  pointer  variable  either  points  to 
the  beginning  of  a  block  of  computer  store  allocated  to  a  node,  or  it  is  nil.  A 
pointer  has  the  value  nil  when  it  points  to  “nothing",  i.e.  it  does  not  point  to  the 
memory  location  of  a  node.  For  example,  the  pointer  “parent”  would  point  to 
the  beginning  of  the  information  for  the  parent  node  of  x tn.  (The  parent  pointer 
of  the  root  would  be  nil.)  The  pointer  “left. child"  would  point  to  its  left  child, 
etc.  The  numeric  variables  contain  such  information  as  the  value  t/”,  the  value 
of  the  integral  (3.1)  on  the  interval  I,  etc.  Figure  2  contains  the  description  of  a 
node 

The  function  /  is  an  auxiliary  function  introduced  to  calculate  the  integral 
(3. 1)  For  all  feR /'(£)=  |/"(£)|.  it  follows  that 


/  U&I  +  \f'(un)\(u?)zdx 


UTn  ~Ur 


(/(OiVi)-/W)). 


(6.1) 


since  V ^  is  0  for  all  xe(x<*,xi,*+i).  For  each  minimal  interval,  use  (6.1)  to  calcu¬ 
late  right. integral  and  left. integral;  fur  each  interior  node,  left. integral  and 
right. integral  are  the  values  of  this. integral  for  the  left  and  right  children  of  the 
node.  For  any  node,  calculate  this. integral  by  adding  the  difference  of  the  left 
and  right  derivatives  of  u  at  the  meshpoint  x?  to  the  sum  of  left  in' egral  and 
right. integral. 


In  each  node  we  have  included  pointers  and  numerical  variables  whose 
values  can  be  caleulatcd  from  already  available  information.  For  example,  in 
each  leaf  node  x"  wc  save  A(*  and  These  values  arc  needed  often  in  the 
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Node: 

x  j  *in  j 
u  jl//*} 

isoleaf  j  true  if  this  node  is  a  leaf  | 
Left. child,  right.child  f  pointers  to 


and  — - - if  this  node  is  not  a 

(C 


parent  {  pointer  to  parent  of  this  node  | 

isoleftchild  {  true  if  this  node  is  a  left  child  } 

depth  |  distance,  in  nodes,  from  this  node  to  the  root  node  | 

left. neighbor.  right.ncighbor  |  pointers  to  zf_,  and  zft.,  ) 

left. boundary,  right. boundary  {  pointers  to  z(  and  zr  ) 

left. sibling.  right.sibling  j  pointers  to  the  nearest  nodes  at  the  same  depth 
as  zf.  to  the  left  and  right,  respectively,  of  zj1  j 
this. integral,  left. integral,  right. integral  {  the  integral  (3.1)  over  (zt,zr), 
(z^z"),  and  (z*\zr),  respectively  | 

,  u&\ -U?  U?-ur-i  , 

'  hi"  1  hj* 

,  t/?M  -  u?  , 

left.ux,  right. ux  |  if  this  node  is  a  leaf,  - — - and  — — - j 

hi  _  hi+1 

fluxplus.  fluxminus,  fluxbar  j  /  v( f/"),  /  “( f/").  /  ( £/")  j 

delta toverh,  hinverse,  criticalvaluc  j  A t/ (zJ*-Zj),  1/  (z^-zt),  c/  (zr-z{)  j 


Figure  2.  Definition  of  a  node  corresponding  to  ztn  in  the  interval  (z,  ,zr). 


scheme,  and  change  only  when  the  mesh  is  changed.  Our  experience  with  the 
scheme  has  shown  that,  on  average,  less  than  one  node  is  added  to  or  sub¬ 
tracted  trom  the  mesh  at  each  time  step.  Thus,  using  the  saved  values  reduces 
execution  time  significantly.  Because  a  new  mesh  is  constructed  at  each  time 
step,  we  also  require  that  the  structural  information  necessary  to  derive  a  new 
mesh  be  readily  available  at  each  node.  Again,  extra  storage  reduces  execution 
time.  Since,  for  many  problems,  many  fewer  nodes  are  necessary  to  obtain  a 
satisfactory  error  with  this  method  than  with  standard  first  order  methods,  the 
extra  storage  requirements  are  not  doomed  critical. 

The  special  nodes  a  and  b  are  the  loft  and  right  boundaries  of  [ a. ,t»  J,  the 
root  node.  In  addition,  supplemental  meshpoints  arc  added  lo  Llie  left  and  to 


-28- 


the  right  of  the  interval  [a,6 ]  so  that  the  pointers  Left. sibling  and  right.sibling 
will  not  be  nil  for  any  node  in  the  tree  headed  by[a,6  ].  The  complete  data  struc¬ 
ture  is  shown  in  Figure  3.  The  lines  in  Figure  3  illustrate  only  the  children  and 
sibling  relationships  of  the  nodes. 


M  I*  M  + 


H-H — I — *- 

U*U) 


Figure  3.  Sibling  and  child  relationships  for  the  mesh  points. 


Figure  4  presents  the  procedure  criteria  that  calculates  the  nodul  values 
needed  to  build  the  new  mesh.  On  each  invocation,  criteria  calculates  the  nodal 
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parameters  for  all  nodes  in  the  subtree  headed  by  p .  Indirection  is  indicated  by 
a  caret  ('*).  The  with  statement  means  that  all  unqualified  nodal  variables  belong 
implicitly  to  the  node  to  which  the  pointer  “p"  points.  For  example, 
“right. neighbor-'. left. ux"  refers  to  p’s  right.neighbor's  left.ux.  Each  time  step, 
the  nodal  values  for  a  and  6  are  calculated  first,  and  then  criteria(root)  is 
called.  The  following  lemma  outlines  a  proof  that  the  procedure  is  correct. 


procedure  criteria(p:  nodepointer); 
begin 

with  pA  do  begin 

fluxbar  :=  fbar(u); 
if  isaleaf  then  begin 

left.ux  :=  (u-  left.neighbor-'.u)  *  hinverse; 

left. integral  :=  abs(left.ux  *  (left. neighbor*. fluxbar  -  fluxbar)); 
right. ux  :=  (right. neighbor-'. u  -  u)  •  hinverse; 
right. integrals  abs(right.ux  •  (right. neighbor*. fluxbar  -  fluxbar)); 
uxx  :=  abs(left.ux  -  right. ux); 

this. integral  :=  left. integral  +  right. integral  +  uxx 
end  else  begin  j  p  does  not  point  to  a  leaf  } 
criteria(left.child); 
c  rite  ria(  right .  child) ; 

right,  integral :  =  right . child''  this . integral; 
left. integral  :=  left. child-'. this. integral; 

uxx  :=  abs(right. neighbor-'. left.ux  -  left. neighbor''. right.ux); 

this,  integral  :=  left  integral  +  right. integral  +  uxx 
end  |  if  isaleaf  then  ...  else  ...  j 
end  |  with  p-'  do  ...  } 
end;  )  criteria  | 


Figure  4.  Procedure  for  calculating  the  integral  (3.1). 


LEMMA  6.1.  Criteria(p)  calculates  the  values  of  right. integral,  left. integral, 
uxx,  and  this  .integral  for  all  nodes  in  the  subtree  headed  by  p.  Furthermore, 
criteria(p)  calculates  left.ux  and  right  .ux  for  all  leaves  in  the  subtree  headed  by 
P 

Proof.  First  one  shows  that  /  (Up)  is  calculated  for  the  left  and  right  boun¬ 
dary  points  of  p  before  criteria(p)  is  called.  This  >s  easily  proved  by  noting  that 
it  is  true  initially  for  the  root  node,  and  that  if  it  is  true  for  p’s  parent,  then  it  is 
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true  for  p.  The  proof  of  the  lemma  then  follows  by  Induction  on  the  height  of  the 
subtree  headed  by  p. 

/// 

For  each  interval  /,  the  integral  (3.1)  over  3/  will  be  used  to  test  for  two 
things: 

1  If  I  is  a  leaf  and  7>e,  to  see  whether  left  and  right  children  of  /  should  be 
added  to  the  set  of  meshpoints,  according  to  the  criteria  of  Step  3  of  Algo¬ 
rithm  M. 

2  If  /  is  not  a  leaf,  to  see  whether  the  subtrees  headed  by  the  left  and  right 
children  of  /  are  still  needed  in  the  mesh. 

These  tests  are  done  on  separate  passes  of  the  tree.  In  the  first  pass,  points  are 
added  to  |x”|  to  obtain  iz"lul2in*1i-  In  the  second  pass,  points  are  removed 
from  the  union  to  leave  only  jxf1*1}.  Since  the  set  of  points  that  are  tested  on 
each  pass  are  for  the  most  part  disjoint  (only  nodes  with  newly  added  children 
are  tested  twice),  this  two  pass  algorithm  does  not  add  too  much  to  the  arith¬ 
metic  complexity  of  the  scheme.  A  one  pass  algorithm  could  surely  be  devised. 

An  argument  similar  to  to  that  in  Lemma  3.1  shows  that,  in  any  mesh  con¬ 
structed  using  Algorithm  M,  a  node's  parent  will  always  have  adjacent  left  and 
right  siblings.  If  jxfl  UW*1!  is  formed  from  jxfj  by  adding  meshpoints  in  order 
of  increasing  depth,  then  this  property  also  holds  for  all  meshes  intermediate  to 
|x/*|  and  |x{*lulxP  +  1l-  Thus,  we  add  nodes  to  the  mesh  in  a  depth  first  ordering, 
and  calculate  (3.1)  for  an  interval  that  is  a  left  child,  for  example,  by  adding 
p~.  parent-',  this,  integral,  p*'. parent*'. left. sibling*'. right. integral,  and 

p~. loft. boundary*'. uxx.  Figure  5  illustrates  the  calculation.  A  similar  expression 
holds  for  nodes  that  are  right  children. 

The  above  calculation  is  the  reason  that  each  admissible  interval  has  a 

’  t 


Figure  S.  The  calculation  of  the  integral  (3.1). 


meshpoint  at  its  center,  and  that  the  integral  (3.1)  is  calculated  over  the  left 
and  right  half  of  each  admissible  interval. 

The  second  pass  of  the  algorithm  checks  each  node  with  children  to  see  if 
the  subtrees  headed  by  that  node's  children  are  needed  to  approximate  the 
function  UJ1*'1  well.  If  not,  they  subtrees  are  removed.  The  nodes  are  checked 
for  subtree  removal  in  a  depth  first  order. 

Note  that  our  algorithm  removes  unnecessary  subtrees  at  once,  instead  of 
removing  one  meshpoint  at  a  time.  In  practice,  however,  we  have  not  not 
observed  the  removal  of  any  but  the  trivial  subtrees  consisting  of  only  one  node. 

We  use  a  simple  and  efficient  memory  management  scheme.  When  a  node  is 
no  longer  needed  and  is  removed  from  the  tree,  it  is  appended  to  a  free  list  of 
nodes.  When  a  node  is  needed  for  addition  to  the  tree,  it  is  taken  from  the  free 
list  if  the  list  is  not  empty.  If  the  free  li3t  is  empty,  a  call  to  the  operating  sys¬ 
tem  allocates  enough  storage  for  the  now  node. 

The  complexity  of  the  scheme  can  be  measured  by  counting  the  number  of 
lloatiug  point  operations  that  are  done,  on  averugc,  for  each  meshpoint  xtn  at 
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each  time  t* .  More  operations  are  performed  for  leaf  nodes  than  for  interior 
nodes.  We  assume  that  an  evaluation  of  /  \  /“,  or  /  requires  two  additions  (or 
subtractions),  two  multiplications,  and  two  comparisons.  This  would  be  the  case, 
for  example,  if  these  functions  were  defined  as  piecewise  quadratic  functions 
over  four  intervals.  Our  assumptions  yield 

Complexity  per  meshpoint  =  17  Additions  +  9.5  Multiplications  A  7  Comparisons. 

This  may  be  compared  with  a  complexity  estimate  of  8  Additions,  4  Multiplica¬ 
tions,  and  4  Comparisons  for  a  careful  implementation  of  a  uniform  mesh  algo¬ 
rithm  with  the  same  spatial  difference  operator  and  Af  -h.  If  an  arbitrary  vari¬ 
able  spaced  mesh  is  chosen  every  timestep,  the  standard  algorithm  will  require 
10  Additions,  5  Multiplications,  2  Divisions,  and  4  Comparisons  per  meshpoint. 
This  does  not  include  whatever  calculations  are  necessary  to  choose  the  mesh. 
Thus,  the  special  placement  of  the  meshpolnts  in  our  algorithm  no  more  than 
doubles  the  work  per  meshpoint. 

Nonarithmetic  operations  must  be  included  in  any  complexity  estimate  of 
these  algorithms.  It  is  more  difficult  to  quantify  this  nonarithmetic  complexity, 
what  may  be  considered  "overhead."  We  give  empirical  results  in  the  next  sec¬ 
tion  that  show  that  the  nonarithmetic  overhead  is  the  same  for  the  fixed  and 
adaptive  mesh  algorithms. 

A  more  serious  difficulty  in  comparing  the  efficiency  of  these  algorithms  is 
that  the  fixed  mesh  algorithm  converges  at  different  rates  for  differing  fluxes  / . 
Wc  will  find  in  the  next  section  that  for  the  problems  for  which  the  fixed  mesh 
algorithm  performs  relatively  well,  the  new  algorithm  compares  poorly.  The 
problems  for  which  the  fixed  mesh  algorithm  performs  poorly,  however,  are 
solved  with  surprising  success  by  our  algorithm. 
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7.  Computational  Results. 

The  Pascal  Implementation,  of  our  algorithm  was  run  on  two  Digital  Equip¬ 
ment  Corporation  VAX  11/760  computers  with  floating  point  accelerators  under 
the  VMS  2.5  and  Berkeley  Unix  4.1  BSD  operating  systems.  Both  programs  use 
double  precision  (64  bit)  floating  point  numbers.  We  only  needed  to  change  real 
constants  from  double  precision  to  single  precision  notation  to  move  the  pro¬ 
gram  from  VMS  to  Unix;  no  other  changes  were  needed. 

We  chose  computational  examples  to  highlight  the  strengths  of  our  method 
as  well  as  point  out  directions  for  improvements.  The  biggest  omission  was  of 
problems  with  smooth  solutions.  (We  did  not  implement  the  algorithm  for 
smooth  solutions  given  in  section  5.)  It  is  easily  shown  that  for  monotone  uni¬ 
form  grid  methods,  the  work  expended  to  achieve  an  accuracy  of  5  at  time  T  is 
proportional  to  8~z.  When  given  a  problem  with  a  smooth  solution,  the  present 
implementation  of  our  method  will  achieve  the  same  accuracy,  but  will  take 
time  at  ^ps  that  are  much  smaller  than  the  average  mesh  spacing.  In  fact,  the 
work  for  our  method  will  be  proportional  to  5~2.  This  shortcoming  will  be  taken 
up  in  a  broader  context  later.  We  have  compared  the  methods  for  problems  with 
shocks,  contact  discontinuities,  and  expansion  waves  following  shocks.  These 
comparisons  lead  us  to  believe  that  our  method  is  superior  to  fixed  mesh  mono¬ 
tone  methods  when  the  discontinuities  in  the  solutions  of  (C)  arc  smeared 
across  more  than  a  fixed  number  of  mesh  intervals  by  the  monotone  schemes. 

In  these  examples,  our  method  is  compared  with  a  finite  difference  scheme 
with  a  fixed,  uniform,  spatial  grid.  This  scheme's  difference  operator  is 

vr'-u?  _ 

to  h  1 1 

Always  At  =h  to  avoid  any  divisions  or  multiplications  in  this  part  of  iho  algo¬ 
rithm.  A  careful  implementation  of  the  comparison  scheme  uses  only  two  func- 


tion  evaluations,  three  subtractions,  and  one  addition  per  meshpoint  per 


timestep. 

For  all  examples  but  the  sixth.  f~—0,  so  that  each  difference  operator 
reduces  to  the  standard  upwind-difference  scheme.  The  initial  values  £/t°  were 
chosen  to  interpolate  the  initial  data  u0  at  the  points  ih.  The  solution  of  the 
Axed  mesh  method  wa3  interpreted  as  a  piecewise  linear  function  taking  on  the 
values  U?  at  the  points  (ih.nbt ).  The  difference  in  ]  between  the  approxi¬ 

mate  solution  and  the  piecewise  linear  interpolant  of  the  true  solution  ■u(x.f) 
was  recorded  as  the  error  of  the  both  schemes. 

Table  1  presents  the  fluxes  /  and  the  initial  and  final  values  of  u  for  the 
computational  experiments.  Table  2  contains  the  errors  and  execution  times, 
respectively,  corresponding  to  various  values  of  A*.  For  the  purposes  of  com¬ 
paring  the  methods,  the  parameter  0.  a  measure  of  the  average  mesh  spacing 
where  the  solution  is  smooth,  is  introduced.  The  parameter  6  has  the  value  V5F 
for  the  adaptive  mesh  method,  and  A*  for  the  fixed  mesh  method.  Table  3  tabu¬ 
lates  an  expression  of  the  error  as  error  =  Ctvnea  and  error  =  C0“  for  each  test. 


TABLE  3 
Error  Decay  Rates 


Test 

Adaotive  Mesh 

Fixed  Mesh 

1 

e  =  .027*“  888 

e  =  .5088' 905 

e  =  .033*' 493 

e  = 

472  6  978 

2 

e  =  .  121  £“388 

e  =  .4190'918 

e  =  .057*"  487 

e  = 

7706  901 

3 

e  =  .377  *~  418 

e  =  1.064 6 1800 

e  =  .506*"  490  * 

e  - 

246  6  995  • 

4 

e  =  .179*"310 

e  =  .6046  9,8 

e  =  .183 1~213 

e  = 

64  3  6  457 

5 

e  =  .002*--373 

e  =  .330  6  909 

e  =  .O01f"371 

e  = 

593  6  738 

6 

e  =  .091*"  381 

e  =  .397  6*  034 

e  =  .007*"  376 

e  = 

6520  749 

*  See  text. 

The  CPU  times  in  Table  2  are  the  "user"  times  reported  by  the  Unix  operat¬ 
ing  system.  The  “system”  time,  which  measures  the  time  spent  by  the  operat¬ 
ing  system  to  service  the  program,  was  not  included.  The  system  tunc  varied 
greatly,  depending  on  the  system  usage,  so  it  was  not  deemed  a  reliable  meas¬ 
ure  of  the  methods'  resource  needs.  The  user  times  corresponded  well  with  the 
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TABLE  1 

Flux.  Initial  and  Final  Values  of  u  (x  .0 


f{u  L 


Unix) 


l/4(ti  +  u8) 

1/  4(n  +  u8) 

1/  2u 

1/  2u 

1  —  (1  — ii)8  if  u  a  1/2 
u 8  if  it  5  1/2 


l-(l-u)8  if « i  1/ 2 
-2u3  +  3uz  — u/  2  if  it  ^1/2 


1  if  *  <  1 
[0  if  x  >  1 

jar-1  if  l<x  £2 
3-ar  if  2<x<3 
0  otherwise 

u  (x +2,4) 

1  if  x  £  1 
[0  if  x  >  1 

1  if  x  £  1 
0  otherwise 


1  if  x  <  1 
|0  otherwise 


1  if  x  <  3 
(0  if  x  >3 

[(x  -2)/3  if  2<x  <2  +  VB 
0  otherwise 

[2 - 4  | x  - 3 1  if  |x  -3|  <  1/2 
0  otherwise 

1  if  x  <  3 
0  if  x  >  3 

1  if  x  <  1 

1  —  (x  —  l)/ 0  if  1«=xs£9-4V5 
0  otherwise 

1  if  x  £  1 

1  —  (x  —  1)/ B  if  l<x<9- 4>/§ 
0  otherwise 


CPU  times  measured  on  the  VMS  operating  system  on  a  lightly  loaded  machine. 
The  error  decay  rates  in  Table  3  are  based  on  a  log-log  least  squares  fitting  of 
the  data.  They  were  calculated  from  the  data  in  Table  2,  but  with  the  full  accu¬ 
racy  of  the  data,  which  is  truncated  in  Table  2. 


Test  1  presents  our  method  in  the  best  light.  This  is  a  Riemann  problem 
with  a  strong  nonlinearity  keeping  the  shock  width  to  within  a  few  of  the  small 
(O(Af))  mesh  intervals.  At  the  same  time,  the  solution  is  constant  outside  the 
shock  region,  so  only  0{- log(Af))  meshpoints  are  necessary  to  accurately 
represent  the  solution.  The  result  is  an  error  of  U(Af)  with  a  computational 
complexity  of  0(-Af  log(Af))- 

Teat  2  is  perhaps  the  worst  comparison  between  the  fixed  mesh  and  adap¬ 
tive  mesh  methods.  The  adaptive  method  gives  no  more  than  accuracy  in 
the  regions  where  the  solution  is  smooth.  To  do  so,  it  uses  U(e-l/z)  meshpoints. 
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TABLE  2 

Computational  Results 


5~l 

r  Test  1. 

Test  2. 

Adaptive  mesh 

Fixed  mesh 

Adaptive  mesh 

Fixed  mesh 

2 

1.252-1 

0.12 

2.250-1 

0.02 

2.468-1 

0.13 

3.419-1 

0.02 

4 

3.237-2 

0.00 

1.252-1 

0.07 

1.174-1 

1.18 

2.079-1 

0.07 

0 

e.  093-3 

5.25 

6.456-2 

0.27 

5.449-2 

1.119-1 

0.20 

16 

2.023-3 

26.33 

3.237-2 

1.07 

3.169-2 

58.85 

6.276-2 

1.07 

32 

5.050-4 

132.65 

1.618-2 

4.37 

1.771-2 

383.15 

2.052-2 

4.48 

64 

1.264-4 

625.60 

8.093-3 

17.65 

9.547-3 

2560.93 

1.300-2 

17.22 

mm 

3.161-5 

Em 

4.046-3 

71.25 

5.129-3 

nHrEErl 

6.702-3 

69.95 

5-1 

Test  3. 

Test  4. 

Fixed  mesh 

Adaptive  mesh 

Fixed  mesh 

2 

0.727-1 

0.12 

1.116-0 

0.02 

3.611-1 

0.10 

4.560-1 

0.02 

4 

4.010-1 

1.22 

8.627-1 

0.05 

1.904-1 

0.08 

3.425-1 

0.07 

0 

1.374-1 

10.02 

6.139-1 

0.25 

1.022-1 

6.45 

2.530-1 

0.23 

16 

5.532-2 

73.72 

3.954-1 

1.00 

5.429-2 

42.57 

1.846-1 

1.00 

32 

2.747-2 

4.02 

2.942-2 

299.60 

1.335-1 

3.98 

64 

1.340-2 

3160.97 

1.229-1 

16.07 

2235.78 

9.592-2 

16.00 

l?fl 

0.091-3 

20316.55 

6.243-2 

64.45 

17124.27  , 

6.860-2 

64.40 

5-1 

Test  5. 

Test  6. 

Adaptive  mesh 

Fixed  mesh 

Adaptive  mesh 

Fixed  mesh 

2 

1.000-1 

0.13 

3.924-1 

0.02 

2.010-1 

0.12 

4.234-1 

4 

0.160-2 

0.93 

1.795-1 

0.07 

9.276-2 

0.92 

1.969-1 

0.07 

0 

4.133-2 

6.92 

1.393-1 

0.27 

4.633-2 

6.45 

1.477-1 

0.30 

16 

2.063-2 

38.90 

7.100-2 

1.12 

2.174-2 

41.32 

7.532-2 

1.10 

32 

1.042-2 

238.97 

5.131-2 

4.45 

1.069-2 

253.40 

5.345-2 

4.68 

64 

5.390-3 

1472.12 

2.875-2 

17.63 

5.440-3 

1594.37 

2.903-2 

120 

2.710-3 

9701.87 

1.590-2 

71.07 

2.716-3 

10267.05 

1.644-2 

75.70 

Coupled  with  a  time  step  ol  0(c),  we  achieve  a  complexity  bound  of  0(e  3/2),  a 


decidedly  inferior  result  when  compared  with  the  fixed  mesh.  This  predicament 
strongly  suggests  a  method  where  the  local  timesteps  are  proportional  to  the 
local  meshsize.  Such  methods  have  been  used  by  Oliger  [  10],  Bolstud  [2],  Berger 
[1 J.  and  others  (see  references  in  [  13j).  With  such  a  scheme,  the  error  will  still 
be  0(e,/z),  but  the  complexity  of  the  scheme  will  be  reduced  to  0(e).  This  is  the 
same  relationship  between  error  size  and  complexity  that  applies  for  the  fixed 
mesh  method,  which  does  surprisingly  well  for  this  problem. 

When  applied  to  the  third  test  problem,  the  fixed  mesh  method  be  haven  in  a 


diflerenL  way  when  the  meshsize  is  small  than  when  it  is  large.  This  is  because  oi 
the  large  second  derivatives  in  the  solution.  The  fixed  mesh  method  '-'iliibits  an 
error  of  order  Af",  with  a  decidedly  less  than  one,  when  the  meshsize  is  large 
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(greater  than  1/  32).  When  the  mesh  size  is  small,  however,  the  mesh  can  ade¬ 
quately  refine  the  large  gradients  and  the  error  is  of  order  A t .  Because  we  are 
mainly  interested  in  exhibiting  “asymptotic"  error  rates,  we  have  computed  the 
error  for  the  fixed  mesh  method  for  At  as  small  as  1/  1024,  and  our  error  decay 
rates  in  Table  3  are  derived  from  the  smaller  timesteps. 

Test  4  was  the  problem  that  motivated  the  design  of  the  adaptive  method. 
It  corresponds  to  a  contact  discontinuity  in  gas  dynamics.  The  discontinuity  is 
smeared  over  0(Af  _l/z)  intervals,  and  the  fixed  mesh  method  has  an  accuracy  of 
order  A/#  The  time  complexity  of  the  fixed  mesh  method  is  0(Af-z),  so  the 
error  of  the  method  is  of  order  t~1/4.  The  adaptive  method  puts  0(Af-I/z)  inter¬ 
vals  of  size  Af  near  the  discontinuity,  and  the  same  number  of  size  At 1/2  away 
from  the  it.  Thus,  the  complexity  of  this  scheme  is  of  order  A f_3/z,  and  the 
error  is  of  order  t~l/a.  In  other  words,  if  one  wants  an  error  of  order  e,  it  takes 
order  c~*  work  for  the  fixed  method,  and  order  e-a  work  for  the  adaptive 
method. 

The  solutions  of  Tests  5  and  6  exhibit  behavior  similar  to  the  solutions  of  the 
Buckley-Leverett  equation  used  as  a  model  in  petroleum  reservoir  simulation 
(see  Douglas  and  Wheeler  [8]).  These  solutions  consist  of  a  shock  followed  by  a 
smooth  expansion  wave.  The  shock  is  not  as  strong  as  in  the  first  example.  The 
fixed  mesh  solution  reflects  this,  as  the  shock  is  smeared  over  several  mesh 
intervals,  and  an  accuracy  of  order  A t3/41  instead  of  At,  is  observed.  The  tests 
indicate  that  the  relationship  between  complexity  and  error  is  the  same  for  the 
adaptive  method  and  the  fixed  mesh  method.  Changing  the  adaptive  method  to 
use  locally  varying  timesteps  as  well  as  mesh  spacing  would  deer  5e  the  CPU 
time  greatly  for  these  problems. 


The  results  of  Tests  1  through  6  suggest  that  our  new  method  is  effective 
when  discontinuities  in  the  solutions  are  smeared  across  many  mesh  intervals 
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by  the  fixed  mesh  method.  This  happens  in  particular  for  problems  that  have 
discontinuous  solutions  but  weak  nonlinear  fluxes. 

An  alternate  approach  to  developing  an  adaptive  mesh  method  is  presented 
in  [6].  There,  the  motivation  is  to  use  an  implicit  scheme  with  a  large  timestep 
and  to  refine  the  mesh  near  roughness  in  the  solution.  Unfortunately,  this  stra¬ 
tegy  does  not  asymptotically  decrease  the  error.  Heuristicaily,  this  is  because 
one  cannot  drag  a  shock  or  discontinuity  across  many  meshpoints  in  one 
timestep  without  smearing  the  discontinuity.  This  effect  is  illustrated  in  Table  4 
for  two  problems.  Each  problem  has  the  solution  u(x)=  1/2  (l-sgn(x-f  /2)). 
The  fluxes  are  f(u)  =  l/2u  and  /  (u)  =  1/ 4  {u  +  u2), respectively.  Each  prob¬ 
lem  was  run  with  h,  the  mesh  spacing,  equal  to  Af  and  Af8.  Thus,  our  choice  is 
not  to  adapt  the  mesh,  but  to  use  a  uniformly  fine  mesh.  By  using  a  fine  mesh 
everywhere,  our  results  do  not  depend  on  any  particular  mesh  refinement  algo¬ 
rithm.  For  the  first  problem  the  error  is  of  order  Af 1/2  for  both  choices  of  the 
mesh.  The  error  for  the  second  method  is  of  order  At.  Thus,  it  appears  that 
spatial  mesh  refinement,  without  a  corresponding  refinement  of  the  temporal 
increments,  is  not  effective  in  solving  these  problems. 


TABLE  4 

Errors  Using  Implicit  Method 


At 

H~  Test  1. 

1  h=Mz  h  =  At 

Test  2. 

h=  At2  h  =  Af 

0 500000 

[  7  »550c-01 

9.5360c-0l 

5. I3l2o-01 

tv  70-'  Be-01 

0  250000 

4.0564e-01 

G.U204C-01 

2.3904e-01 

4. 1525e-0l 

0.125000 

3.1453c-01 

4.05G4c-0l 

1.0772e-01 

2.3833e-01 

o  oir ’!>oo 

21 i20c-01 

3.4  ^  *e-01 

4.9120c-02 

•  .27** *  e-P* 

o  03’  :t  ->o 

1 .4!>29c-01 

2.4393o-01 

2. 3225e-02 

G.4799e-i>2 

0  015625 

l  1  0124C-01 

1.720  lc-01 

1 . 1274e-02 

3.2443e-02 

In  the  previous  section,  we  showed  that  the  arithmetic  complexity  of  the 
adaptive  .  tlgorithm  was  no  more  than  twice  that  of  the  fixed  point  algo- 

n tlun  per  mcsiipoiiit  per  timestep.  The  algorithms  diner  _  realiy  u.  their  imple- 
nteiii  itioiij,  however  The  fixed  mesh  algorithm  is  almost  trivial  to  implement, 
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while  the  adaptive  method  uses  sophisticated  data  structures  and  pointer  mani¬ 
pulation.  We  therefore  set  out  to  quantify  the  amount  of  nonarithmetic  over¬ 
head  in  each  method. 

We  proposed  to  measure  the  proportion  of  the  CPU  time  spent  in  arithmetic 
computations  as  compared  with  nonarithmetic  computations  for  both  imple¬ 
mentations.  We  had  available  two  VAX’s  running  identical  software,  one  of  which 
did  not  have  a  floating  point  accelerator  (FPA).  A  simple  program  consisting 
mainly  of  an  equal  number  of  nontrivial  memory  to  register  floating  point  addi¬ 
tions  and  multiplications  was  run  and  timed  on  both  machines.  A  speedup  of  a 
factor  of  five  was  observed  on  the  machine  with  the  FPA.  (The  effects  of  cache 
memory  hits  was  not  considered.)  We  then  compared  the  CPU  times  for  the 
same  implementations  of  the  two  algorithms  on  both  machines.  The  ratios  of 
the  CPU  times  are  given  in  Table  5.  Assuming  that  only  floating  point  operations 
were  speeded  up  by  the  FPA  (a  faulty  assumption,  for  some  integer  arithmetic  is 
also  faster),  we  calculated  the  fraction  of  time  spent  in  nonarithmetic  opera¬ 
tions  by  both  programs  on  both  machines.  For  Test  6,  with  nontrivial  /*,  f~, 
and  / ,  around  90%  of  the  time  was  spent  on  nonarithmetic  operations  on  the 
machine  with  the  FPA  for  both  problems.  The  figure  is  similar  for  Test  4,  with 
trivial  /  \  /  ”,  and  / . 

To  discover  the  effects  of  the  Pascal  compiler,  we  implemented  the  fixed 
mesh  algorithm  in  FORTRAN  and  ran  this  program  on  the  machine  with  the  FPA. 
For  sonic  reason,  the  FORTRAN  compiler  generated  much  superior  code  for 
index  calculations.  Function  calls  were  slightly  cheaper  because  the  FORTRAN 
program,  not  being  block  structured,  did  not  maintain  a  "display”  of  the 
currently  accessible  data  areas.  The  computation  time  of  the  FORTRAN  program 
was  t>d%  of  its  Pascal  counterpart.  This  still  leaves  us  with  an  overhead  figure  of 
about  07%.  These  tests  indicate  that  the  overhead  is  similar,  and  large,  for  each 


algorithm. 


TABLE  5 

NONAniTHMETIC  OVERHEAD 


Adaptive  mesh 

Fixed  mesh 

Test 

CPU  time 

Overhead 

CPU  time 

Overhead 

ratio 

FPA  No  FPA 

ratio 

FPA 

No  FPA 

4 

0.728 

91%  66% 

0.678 

88% 

60% 

6 

0.699 

89%  63% 

0.715 

90% 

64% 

B.  REFERENCES. 

1.  M.  Berger,  "Adaptive  mesh  refinement  for  hyperbolic  partial  differential 
equations."  Stanford  Computer  Science  Report  STAN-CS-82-924  (disserta¬ 
tion). 

2.  J.  H.  Bolstad,  "An  adaptive  finite  difference  method  for  hyperbolic  systems 
in  one  space  dimension,"  Lawrence  Berkeley  Lab.  LBL-132B7  (STAN-CS-B2- 
899)  (dissertation). 

3.  C.  de  Boor.  “Good  approximation  by  splines  with  variable  knots,"  in  Spline 
functions  and  approximation  theory,  A  Meir  and  A.  Sharma  ed  ,  J3NM  v.  21, 
Birkhauser  Verlag,  1973,  pp.  57-72. 

4.  C.  de  Boor,  "Good  approximation  by  splines  with  variable  knots  II,”  in  Lec¬ 
ture  Notes  in  Mathematics  363,  Springer  Verlag,  1974,  pp.  12-20. 

5.  M.  G.  Crandall  Sc  A.  Majda,  "Monotone  difference  approximations  for  scalar 
conservation  laws,"  Math.  Comp  v.  34,  1900,  pp.  1-21. 

6.  S.  F.  Davis  Sc  J.  E.  Flaherty,  "An  adaptive  finite  element  method  for  initial- 
value  problems  for  partial  differential  equations,"  SIAM  J.  Sci.  Slat.  Corn- 
put.,  v.  3,  1902,  pp,  6-20. 

7.  J.  Douglas  Jr.,  "Simulation  of  a  linear  waterflood,”  in  I'Yee  Uoundary  Jh-ob- 
lems,  Proceedings  of  a  seminar  held  in  Pavia,  Sept. -Oct.  1979,  Vol.  II,  Insti- 

i 


-39- 


tuto  Nazionale  dl  Alta  Matematica  "Francesco  Seven."  Roma,  1900. 

6.  J.  Douglas  Jr.  &  M.  F.  Wheeler,  "Implicit,  time-dependent  variable  grid  finite 
difference  methods  for  the  approximation  of  a  linear  waterflood,"  Math. 
Comp.,  v.  40,  1983,  pp.  107-122. 

9.  Todd  Dupont,  "Mesh  modification  for  evolution  equations,"  Math.  Comp.,  v. 
39,  1902,  pp.  85-107. 

10.  B.  Enquist  &  S.  Osher,  "Stable  and  entropy  satisfying  approximations  for 
transonic  flow  calculations.”  Math.  Comp.,  v.  34.  1980,  pp.  45-75. 

11.  D.  B.  Gannon,  "Self  adaptive  methods  for  parabolic  partial  differential  equa¬ 
tions."  U.  of  Illinois  CS  Dept,  report  UIUCDCS-R-80-1020. 

12.  A.  Harten,  J.  M.  Hyman  &  P.  D.  Lax,  "On  finite  difference  approximations  and 
entropy  conditions  for  shocks,"  Comm.  Pure  Appl.  Math.,  v.  29.  1976,  pp. 
297-322. 

13.  G.  W.  Hedstrom  &  G.  H.  Rodrigue,  "Adaptive-grid  methods  for  time- 
dependent  partial  differential  equations,”  in  Multigrid  Methods,  W.  Hack- 
busch  and  U.  Trottenberg,  ed..  Springer  Verlag,  1982,  pp.  474-404. 

1  S.  N.  Kruzkov,  "First  order  quasilincar  equations  with  several  space  vari¬ 
ables,"  Math.  USSR  Sb.,  v.  10,  1970,  pp.  217-243. 

15.  N.  N.  Kuznetsov,  “On  stable  methods  for  solving  non-linear  first  order  par¬ 
tial  differential  equations  in  the  class  of  discontinuous  functions,"  in  Topics 
in  Numerical  Analysis,  John  J.  H.  Miller,  cd.,  Academic  Press,  New  York, 
1977,  pp.  183-197. 

16.  B.  J.  Lucier,  "On  nonlocal  monotone  difference  methods  for  scalar  conserva¬ 
tion  laws,"  to  appear. 

17  K.  Miller,  "Moving  finite  elements,  parts  II,”  SIAM  J.  Numer.  Anal.,  v.  10, 
1981,  pp  1019-1057. 


-40- 


18.  J.  Oliger,  “Approximate  Methods  for  Atmospheric  and  Oceanographic  Circu¬ 
lation  Problems,"  in  lecture  Notes  in  Physics  91,  R.  Glowinski  and  J.  Lions, 
ed..  Springer  Verlag.  1979,  pp.  171-184. 

19.  R.  Sanders,  "On  convergence  of  monotone  finite  difference  schemes  with 
variable  spatial  differencing,”  Math.  Comp.,  v.  40.  1983,  pp.  91-106. 

i 


BJL/jvs 


SECURITY  CLASSIFICATION  OF  THIS  PAOC  fWtMn  S«ifaMn« 


|  REPORT  DOCUMENTATION  PACE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

S.  RECIPIENT'S  CATALOG  NUMGER 

4.  TITLE  (end  Subtitle) 

A  Stable  Adaptive  Numerical  Scheme  for 

Hyperbolic  Conservation  laws 

Summary  Report  -  no  specific 
reporting  period 

«.  PERFORMING  ORO.  REPORT  NUMBER 

7.  AUTMORfA) 

Bradley  J.  Lucier 

i.  Contract  or  grant  number^ 

MCS-7927062,  Mod.  1 
DAAG2§l80-6-0041 

PERFORMING  ORGANIZATION  NAME  AND  ADORESS 

Mathematics  Research  Center,  University  of 

610  Walnut  Street  Wisconsin 

Madison.  Wisconsin  53706 

to.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  A  WORK  UNIT  NUMBERS 

Work  Unit  Number  3  - 
Numerical  Analysis  and 
Scientific  Comoutina 

II.  CONTROLLING  OFFICE  NAME  ANO  ADORESS 

Sea  Item  18  below 

IZ.  REPORT  oate 

May  1983 

IS.  NUMBER  OF  PAGES 

40 

14-  MONITORING  AOENCY  NAME  A  AODRESSflf  dltletmH  bum  CentrolUnd  Oil lee) 

IS.  SECURITY  CLASS,  (o 1  thle  neon) 

UNCLASSIFIED 

IS.  DISTRIBUTION  STATEMENT  (el  Bile  mOm I) 


Approved  for  public  release;  distribution  unlimited. 


17.  DISTRIBUTION  STATEMENT  (•(  AM  obmUeel  fat  Block  so,  II  Ollmil  tw  Jhwn) 


IS.  SUPPLEMENTARY  NOTES 

U.  S.  Army  Research  Office  National  Science  Foundation 

P.  0.  Box  12211  Washington,  DC  20550 

Research  Triangle  Park 
North  Carolina  27709 

It.  KEY  WORDS  (Continue  on  tewetee  oldo  II  nmoorrery  end  Identity  by  Hook  number)  ’ 

Conservation  laws,  finite-difference  schemes,  adaptive  numerical  methods. 


SO.  ABSTRACT  (Continue  on  rerot ee  eld e  II  neeeee ety  mtd  Identity  by  block  number) 

A  new  adaptive  finite-difference  scheme  for  scalar  byperholic  conservation 
laws  is  introduced.  A  key  aspect  of  the  method  is  a  new  automatic  mesh 
selection  algorithm  for  problems  with  shocks.  We  show  that  the  scheme  is  IT- 
stable  in  the  sense  of  Kuznetsov,  and  that  it  generates  convergent  approximations 
for  linear  problems.  Numerical  evidence  is  presented  that  indicates  that  if  an 
error  of  size  e  is  required,  our  scheme  takes  at  most  0(e“^)  operations. 
Standard  monotone  difference  schemes  can  take  up  to  0(c~4)  calculations  for  the 
same  problems. 
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